In [1]:
import importlib
import theano.tensor as T
import sys, os
sys.path.append("/home/bl3/PycharmProjects/GeMpy/")
sys.path.append("/home/bl3/PycharmProjects/pygeomod/pygeomod")
import GeoMig
#import geogrid
#importlib.reload(GeoMig)
importlib.reload(GeoMig)
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import matplotlib
import matplotlib.cm as cmx
from skimage import measure
os.environ['CUDA_LAUNCH_BLOCKING'] = '1'
np.set_printoptions(precision = 2, linewidth= 130, suppress = True)
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
%matplotlib qt4
In [2]:
interpolator = GeoMig.GeoMigSim_pro2(c_o = np.float32(-0.1),range = np.float32(110000)) # Range used in geomodeller
# Setting the extent and resolution of the grid
xmin = 696000
xmax = 747000
ymin = 6863000
ymax = 6950000
zmin = -20000
zmax = 2000
"""
xmin = layers[:,0].min()
xmax = layers[:,0].max()
ymin = layers[:,1].min()
ymax = layers[:,1].max()
zmin = layers[:,2].min()
zmax = layers[:,2].max()
"""
nx = 50
ny = 50
nz = 50
dx = (xmax-xmin)/nx
dy = (ymax-ymin)/ny
dz = (zmax-zmin)/nz
#Genereting a regular grid
interpolator.create_regular_grid_3D(xmin,xmax,ymin,ymax,zmin,zmax,nx,ny,nz)
interpolator.theano_set_3D_nugget()
In [42]:
import pandas as pn
Foliations = pn.read_csv("a_Foliations.csv")
Points = pn.read_csv("a_Points.csv")
Foliations[Foliations["formation"] == "EarlyGranite"];
So there is four series, 3 of one single layer and 1 with 2. Therefore we need 4 potential fields, so lets begin.
In [4]:
# We extract the points of the formations and we make a numpy array float 32 out of it
layers = Points[Points["formation"] == "EarlyGranite"].as_matrix()[:,:-1].astype("float32")
layers
# The code for rest and ref has to change a bit according to the number of layer per formation
rest = layers[1:]
ref = np.tile(layers[0],(np.shape(layers)[0]-1,1))
rest, ref;
In [5]:
dips = Foliations[Foliations["formation"] == "EarlyGranite"].as_matrix()[:,:3].astype("float32")
dips_angles = Foliations[Foliations["formation"] == "EarlyGranite"]["dip"].as_matrix().astype("float32")
azimuths = Foliations[Foliations["formation"] == "EarlyGranite"]["azimuth"].as_matrix().astype("float32")
polarity = Foliations[Foliations["formation"] == "EarlyGranite"]["polarity"].as_matrix().astype("float32")
In [6]:
G_x = np.sin(np.deg2rad(dips_angles)) * np.sin(np.deg2rad(azimuths)) * polarity
G_y = np.sin(np.deg2rad(dips_angles)) * np.cos(np.deg2rad(azimuths)) * polarity
G_z = np.cos(np.deg2rad(dips_angles)) * polarity
G_x, G_y, G_z;
In [7]:
b = interpolator.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0]
arg1,h1 = np.argmin((abs(interpolator.grid - ref[0])).sum(1)), b[np.argmin((abs(interpolator.grid - ref[0])).sum(1))]
arg2, h2 =np.argmin((abs(interpolator.grid - ref[1])).sum(1)), b[np.argmin((abs(interpolator.grid - ref[1])).sum(1))]
In [8]:
xmin = layers[:,0].min()
xmax = layers[:,0].max()
ymin = layers[:,1].min()
ymax = layers[:,1].max()
zmin = layers[:,2].min()
zmax = layers[:,2].max()
nx = 50
ny = 50
nz = 50
dx = (xmax-xmin)/nx
dy = (ymax-ymin)/ny
dz = (zmax-zmin)/nz
#Genereting a regular grid
interpolator.create_regular_grid_3D(xmin,xmax,ymin,ymax,zmin,zmax,nx,ny,nz)
interpolator.theano_set_3D_nugget()
val = interpolator.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0]
G_x, G_y, G_z = interpolator.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[-3:]
sol_early_granite = val.reshape(50,50,50,
order = "C")
arg1,h1 = np.argmin((abs(interpolator.grid - ref[0])).sum(1)), val[np.argmin((abs(interpolator.grid - ref[0])).sum(1))]
arg2, h2 =np.argmin((abs(interpolator.grid - ref[1])).sum(1)), val[np.argmin((abs(interpolator.grid - ref[1])).sum(1))]
In [34]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
cm = plt.get_cmap("jet")
cNorm = matplotlib.colors.Normalize(vmin=sol_early_granite.min(), vmax=sol_early_granite.max())
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)
isolines = np.linspace(h1,sol_early_granite.min()+(sol_early_granite.max()-sol_early_granite.min())*0.01,4)
print(isolines)
for i in isolines[0:10]:
vertices = measure.marching_cubes(sol_early_granite, i, spacing = (dx,dy,dz),
gradient_direction = "descent")[0]
# ax.scatter(vertices[::10,0]+xmin,vertices[::10,1]+ymin,zmin+vertices[::10,2],color=scalarMap.to_rgba(i),alpha = 0.2)
#color=scalarMap.to_rgba(vertices[::10,2])
ax.scatter(layers[:,0],layers[:,1],layers[:,2], s = 50, c = "r" )
ax.scatter(ref[0,0],ref[0,1],ref[0,2], s = 100, c = "g")
ax.quiver3D(dips[:,0],dips[:,1],dips[:,2], G_x*size,G_y*size,G_z*size, pivot = "tail", length = 40,
arrow_length_ratio = 0.01 )
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.ticklabel_format(style='sci', scilimits=(0,0))
In [30]:
In [75]:
# We extract the points of the formations and we make a numpy array float 32 out of it
layer_1 = Points[Points["formation"] == "SimpleMafic2"].as_matrix()[:,:-1].astype("float32")
layer_2 = Points[Points["formation"] == "SimpleBIF"].as_matrix()[:,:-1].astype("float32")
layers = np.asarray([layer_1,layer_2])
# The code for rest and ref has to change a bit according to the number of layer per formation
rest = np.vstack((i[1:] for i in layers))
ref = np.vstack((np.tile(i[0],(np.shape(i)[0]-1,1)) for i in layers))
"""
dips = np.asarray(
[Foliations[Foliations["formation"] == "SimpleMafic2"].as_matrix()[:,:3].astype("float32"),
Foliations[Foliations["formation"] == "SimpleBIF"].as_matrix()[:,:3].astype("float32")]
)
dips_angles = np.asarray(
[Foliations[Foliations["formation"] == "EarlyGranite"]["dip"].as_matrix().astype("float32"),
Foliations[Foliations["formation"] == "EarlyGranite"]["dip"].as_matrix().astype("float32")]
)
azimuths = np.asarray(
[Foliations[Foliations["formation"] == "SimpleMafic2"]["azimuth"].as_matrix().astype("float32"),
Foliations[Foliations["formation"] == "SimpleBIF"]["dip"].as_matrix().astype("float32")]
)
polarity = np.asarray(
[Foliations[Foliations["formation"] == "SimpleMafic2"]["polarity"].as_matrix().astype("float32"),
Foliations[Foliations["formation"] == "SimpleBIF"]["polarity"].as_matrix().astype("float32")]
)
"""
dips = np.asarray(
Foliations[Foliations["formation"] == "SimpleBIF"].as_matrix()[:,:3].astype("float32")
)
dips_angles = np.asarray(
Foliations[Foliations["formation"] == "SimpleBIF"]["dip"].as_matrix().astype("float32")
)
azimuths = np.asarray(
Foliations[Foliations["formation"] == "SimpleBIF"]["azimuth"].as_matrix().astype("float32")
)
polarity = np.asarray(
Foliations[Foliations["formation"] == "SimpleBIF"]["polarity"].as_matrix().astype("float32")
)
ref
Out[75]:
In [90]:
xmin = 696000
xmax = 747000
ymin = 6863000
ymax = 6950000
zmin = -20000
zmax = 2000
"""
xmin = layers[0][:,0].min()
xmax = layers[0][:,0].max()
ymin = layers[0][:,1].min()
ymax = layers[0][:,1].max()
zmin = layers[0][:,2].min()
zmax = layers[0][:,2].max()
"""
nx = 50
ny = 50
nz = 50
dx = (xmax-xmin)/nx
dy = (ymax-ymin)/ny
dz = (zmax-zmin)/nz
#Genereting a regular grid
interpolator.create_regular_grid_3D(xmin,xmax,ymin,ymax,zmin,zmax,nx,ny,nz)
In [94]:
interpolator.theano_set_3D_nugget()
val = interpolator.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0]
G_x, G_y, G_z = interpolator.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[-3:]
print()
sol_early_granite = val.reshape(50,50,50,
order = "C")
arg1,h1 = np.argmin((abs(interpolator.grid - ref[0])).sum(1)), val[np.argmin((abs(interpolator.grid - ref[0])).sum(1))]
arg2, h2 =np.argmin((abs(interpolator.grid - ref[-1])).sum(1)), val[np.argmin((abs(interpolator.grid - ref[-1])).sum(1))]
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
cm = plt.get_cmap("jet")
cNorm = matplotlib.colors.Normalize(vmin=sol_early_granite.min(),
vmax=sol_early_granite.max())
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)
isolines = np.linspace(sol_early_granite.min()+(sol_early_granite.max()-sol_early_granite.min())*0.01,
sol_early_granite.max()-(sol_early_granite.max()-sol_early_granite.min())*0.01,
10)
print(isolines)
for i in isolines[0:10]:
vertices = measure.marching_cubes(sol_early_granite, i, spacing = (dx,dy,dz),
gradient_direction = "descent")[0]
ax.scatter(vertices[::10,0]+xmin,vertices[::10,1]+ymin,zmin+vertices[::10,2],color=scalarMap.to_rgba(i),alpha = 0.2)
#color=scalarMap.to_rgba(vertices[::10,2])
ax.scatter(layers[0][:,0],layers[0][:,1],layers[0][:,2], s = 50, c = "r" )
ax.scatter(layers[1][:,0],layers[1][:,1],layers[1][:,2], s = 50, c = "yellow" )
ax.scatter(ref[0,0],ref[0,1],ref[0,2], s = 100, c = "g")
ax.quiver3D(dips[:,0],dips[:,1],dips[:,2], G_x,G_y,G_z, pivot = "tail", length = 400,
arrow_length_ratio = 0.01 )
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.ticklabel_format(style='sci', scilimits=(0,0))
In [82]:
h1,h2, sol_early_granite.min(), sol_early_granite.max()
Out[82]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [14]:
"""Export model to VTK
Export the geology blocks to VTK for visualisation of the entire 3-D model in an
external VTK viewer, e.g. Paraview.
..Note:: Requires pyevtk, available for free on: https://github.com/firedrakeproject/firedrake/tree/master/python/evtk
**Optional keywords**:
- *vtk_filename* = string : filename of VTK file (default: output_name)
- *data* = np.array : data array to export to VKT (default: entire block model)
"""
vtk_filename = "noddyFunct2"
extent_x = 10
extent_y = 10
extent_z = 10
delx = 0.2
dely = 0.2
delz = 0.2
from pyevtk.hl import gridToVTK
# Coordinates
x = np.arange(0, extent_x + 0.1*delx, delx, dtype='float64')
y = np.arange(0, extent_y + 0.1*dely, dely, dtype='float64')
z = np.arange(0, extent_z + 0.1*delz, delz, dtype='float64')
# self.block = np.swapaxes(self.block, 0, 2)
gridToVTK(vtk_filename, x, y, z, cellData = {"geology" : sol})
In [ ]:
In [ ]:
In [ ]:
In [32]:
%%timeit
sol = interpolator.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0]
In [17]:
interpolator.geoMigueller.profile.summary()
In [ ]:
In [23]:
sys.path.append("/home/bl3/anaconda3/lib/python3.5/site-packages/PyEVTK-1.0.0-py3.5.egg_FILES/pyevtk")
nx = 50
ny = 50
nz = 50
xmin = 1
ymin = 1
zmin = 1
grid = sol
var_name = "Geology"
#from evtk.hl import gridToVTK
import pyevtk
from pyevtk.hl import gridToVTK
# define coordinates
x = np.zeros(nx + 1)
y = np.zeros(ny + 1)
z = np.zeros(nz + 1)
x[1:] = np.cumsum(delx)
y[1:] = np.cumsum(dely)
z[1:] = np.cumsum(delz)
# plot in coordinates
x += xmin
y += ymin
z += zmin
print (len(x), x)
gridToVTK("GeoMigueller", x, y, z,
cellData = {var_name: grid})
Out[23]:
In [23]:
interpolator.theano_set_3D()
In [26]:
%%timeit
sol = interpolator.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0].reshape(20,20,20, order = "C")
In [27]:
interpolator.geoMigueller.profile.summary()
In [19]:
interpolator.theano_set_3D_legacy()
In [21]:
%%timeit
sol = interpolator.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0].reshape(20,20,20, order = "C")
In [22]:
interpolator.geoMigueller.profile.summary()
In [1]:
from theano import function, config, shared, sandbox
import theano.tensor as T
import numpy
import time
vlen = 10 * 30 * 768 # 10 x #cores x # threads per core
iters = 1000
rng = numpy.random.RandomState(22)
x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
f = function([], T.exp(x))
print(f.maker.fgraph.toposort())
t0 = time.time()
for i in range(iters):
r = f()
t1 = time.time()
print("Looping %d times took %f seconds" % (iters, t1 - t0))
print("Result is %s" % (r,))
if numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()]):
print('Used the cpu')
else:
print('Used the gpu')
In [1]:
from theano import function, config, shared, sandbox
import theano.tensor as T
import numpy
import time
vlen = 10 * 30 * 768 # 10 x #cores x # threads per core
iters = 1000
rng = numpy.random.RandomState(22)
x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
f = function([], T.exp(x))
print(f.maker.fgraph.toposort())
t0 = time.time()
for i in range(iters):
r = f()
t1 = time.time()
print("Looping %d times took %f seconds" % (iters, t1 - t0))
print("Result is %s" % (r,))
if numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()]):
print('Used the cpu')
else:
print('Used the gpu')
In [18]:
from theano import function, config, shared, sandbox
import theano.tensor as T
import numpy
import time
vlen = 10 * 30 * 768 # 10 x #cores x # threads per core
iters = 1000
rng = numpy.random.RandomState(22)
x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
f = function([], T.exp(x))
print(f.maker.fgraph.toposort())
t0 = time.time()
for i in range(iters):
r = f()
t1 = time.time()
print("Looping %d times took %f seconds" % (iters, t1 - t0))
print("Result is %s" % (r,))
if numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()]):
print('Used the cpu')
else:
print('Used the gpu')
In [ ]:
In [759]:
np.set_printoptions(precision=2)
test.geoMigueller(dips,dips_angles,rest, ref)[1]
Out[759]:
In [751]:
T.fill_diagonal?
In [758]:
import matplotlib.pyplot as plt
% matplotlib inline
dip_pos_1_v = np.array([np.cos(np.deg2rad(dip_angle_1))*1,
np.sin(np.deg2rad(dip_angle_1))]) + dip_pos_1
dip_pos_2_v = np.array([np.cos(np.deg2rad(dip_angle_2))*1,
np.sin(np.deg2rad(dip_angle_2))]) + dip_pos_2
plt.arrow(dip_pos_1[0],dip_pos_1[1], dip_pos_1_v[0]-dip_pos_1[0],
dip_pos_1_v[1]-dip_pos_1[1], head_width = 0.2)
plt.arrow(dip_pos_2[0],dip_pos_2[1],dip_pos_2_v[0]-dip_pos_2[0],
dip_pos_2_v[1]-dip_pos_2[1], head_width = 0.2)
plt.plot(layer_1[:,0],layer_1[:,1], "o")
plt.plot(layer_2[:,0],layer_2[:,1], "o")
plt.plot(layer_1[:,0],layer_1[:,1], )
plt.plot(layer_2[:,0],layer_2[:,1], )
plt.contour( sol.reshape(50,50) ,30,extent = (0,10,0,10) )
#plt.colorbar()
#plt.xlim(0,10)
#plt.ylim(0,10)
plt.title("GeoBulleter v 0.1")
print (dip_pos_1_v, dip_pos_2_v, layer_1)
In [443]:
n = 10
#a = T.horizontal_stack(T.vertical_stack(T.ones(n),T.zeros(n)), T.vertical_stack(T.zeros(n), T.ones(n)))
a = T.zeros(n)
print (a.eval())
#U_G = T.horizontal_stack(([T.ones(n),T.zeros(n)],[T.zeros(n),T.ones(n)]))
In [6]:
T.stack?ö+aeg
In [ ]:
_squared_euclidean_distances2 = T.sqrt(
(dips ** 2).sum(1).reshape((dips.shape[0], 1)) + (aux_Y ** 2).sum(1).reshape(
(1, aux_Y.shape[0])) - 2 * dips.dot(aux_Y.T))
_squared_euclidean_distances3 = T.sqrt(
(dips ** 2).sum(1).reshape((dips.shape[0], 1)) + (aux_X ** 2).sum(1).reshape(
(1, aux_X.shape[0])) - 2 * dips.dot(aux_X.T))
h3 = T.vertical_stack(
(dips[:, 0] - aux_Y[:, 0].reshape((aux_Y[:, 0].shape[0], 1))).T,
(dips[:, 1] - aux_Y[:, 1].reshape((aux_Y[:, 1].shape[0], 1))).T
)
h4 = T.vertical_stack(
(dips[:, 0] - aux_X[:, 0].reshape((aux_X[:, 0].shape[0], 1))).T,
(dips[:, 1] - aux_X[:, 1].reshape((aux_X[:, 1].shape[0], 1))).T)
r_3 = T.tile(_squared_euclidean_distances2, (2, 1)) # Careful with the number of dimensions
r_4 = T.tile(_squared_euclidean_distances3, (2, 1)) # Careful with the number of dimensions
_ans_d1_3 = (r_3 < self.a) * (
-7 * (self.a - r_3) ** 3 * r_3 * (8 * self.a ** 2 + 9 * self.a * r_3 + 3 * r_3 ** 2) * 1)
/ (4 * self.a ** 7)
_ans_d1_4 = (r_4 < self.a) * (
-7 * (self.a - r_4) ** 3 * r_4 * (8 * self.a ** 2 + 9 * self.a * r_4 + 3 * r_4 ** 2) * 1)
/ (4 * self.a ** 7)
_C_GI = (h3 / r_3 * _ans_d1_3 - h4 / r_4 * _ans_d1_4).T
self._f_CGI = theano.function([dips, aux_X, aux_Y], _C_GI)
In [67]:
import geomodeller_xml_obj
data = geomodeller_xml_obj.GeomodellerClass()
data.load_geomodeller_file('/home/bl3/PycharmProjects/Temp_SandstoneCopy/Temp_SandstoneCopy.xml')
section_elements = data.get_sections() # This reads geomodeller sections. Each section one object
data.get_formation_point_data(section_elements[0])
points_elements = [data.get_formation_point_data(i) for i in section_elements]
print (points_elements[0][0])
points_coord = data.get_point_coordinates(points_elements[0][0])
importlib.reload(geomodeller_xml_obj)
data = geomodeller_xml_obj.GeomodellerClass()
data.load_geomodeller_file('/home/bl3/PycharmProjects/Temp_SandstoneCopy/Temp_SandstoneCopy.xml')
section_dict = data.create_sections_dict()
contact_points_dict = {}
foliation_dict = {}
for i in range(len(section_dict)):
print ("\n\n\n", next (iter (section_dict.values())), "\n")
print ("Elements and their ID \n")
contact_points = data.get_formation_point_data( next (iter (section_dict.values())))
try:
for contact_point in contact_points:
contact_points_dict[contact_point.get("ObservationID")] = contact_point
print (contact_point, contact_point.get("ObservationID"))
except TypeError:
print ("No contact points in the section")
#ObsID = contact_points.get("ObservationID")
foliations = data.get_foliations( next (iter (section_dict.values())))
try:
for foliation in foliations:
# dictionary to access with azimth name
foliation_dict[foliation.get("ObservationID")+"_a"] = foliation
# dictionary to access with dip name
foliation_dict[foliation.get("ObservationID")+"_d"] = foliation
print (foliation, foliation.get("ObservationID"))
except TypeError:
print ("No foliation in the section")
try:
coord_interface = data.get_point_coordinates(contact_points)
except TypeError:
print ("Element does not have iterable objects")
print ("\nDictionaries:\n ", contact_points_dict, "\n", foliation_dict)
print ("\n Contact points", contact_points, "\n", coord_interface, "\n")
print ("foliations" , foliations, "\n")
try:
for i in range(len(foliations)):
print ("azimut:",data.get_foliation_azimuth(foliations[i]))
print ("dip",data.get_foliation_dip(foliations[i]))
print ("coordinates", data.get_foliation_coordinates(foliations[i]))
except TypeError:
print ("No foliation in the section")